Computational Difficulty of Global Variations in the Density Matrix Renormalization Group 
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The density matrix renormalization group (DMRG) approach is arguably the most successful method to nu- 
merically find ground states of quantum spin chains. It amounts to iteratively locally optimizing matrix-product 
states, aiming at better and better approximating the true ground state. To date, both a proof of convergence 
to the globally best approximation and an assessment of its complexity are lacking. Here we establish a result 
on the computational complexity of an approximation with matrix-product states: The surprising result is that 
when one globally optimizes over several sites of local Hamiltonians, avoiding local optima, one encounters in 
the worst case a computationally difficult NP-hard problem (hard even in approximation). The proof exploits a 
novel way of relating it to binary quadratic programming. We discuss intriguing ramifications on the difficulty 
of describing quantum many-body systems. 

PACS numbers: 



"How difficult is it to describe quantum systems in clas- 
sical terms"? This question in its various variants has man- 
ifold implications to several fields of theoretical physics: to 
the context of numerically studying many-body systems of 
condensed-matter physics in their ground-state properties, to 
the question of the superiority of a quantum compared to a 
classical computer, and others. Recently, a renewed interest in 
questions of classically simulating quantum many-body sys- 
tems [1-3] gave rise to a number of new results and simula- 
tion methods, a large number of them motivated or in their 
approach by ideas of quantum information theory [2-8]. The 
arguable workhorse of numerically finding ground states of 
many -body systems, the DMRG method [1-3], was recently 
reassessed and in some ways improved. It seems fair to say 
that the problem of finding ground states of systems with pe- 
riodic boundary conditions or higher-dimensional systems is 
now much better understood than not very long ago [6] . The 
performance of DMRG-type methods has also been quantita- 
tively related to entanglement scaling in ground states: One 
should expect an approximation in terms of matrix-product 
states - as DMRG is generating - to be most faithful, if an 
"entanglement area- theorem" holds, which in turn is typically 
the case in non-critical systems [8, 10]. 
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Now, even if matrix product states form the right set of 
states that well-describe the true ground state properties, and 
one can expect MPS to faithfully represent the ground state 
'How difficult is it then to find the truly best approxi- 
mation to the ground state"? This is the key question of the 
computational complexity of any method to find ground states. 
Quite surprisingly, given the maturity of the field and the sig- 
nificance of such simulations, this question is essentially open. 
In practice, DMRG produces very good results, despite of the 
possibility of getting stuck in local minima in the optimiza- 
tion. However, it is not certifiable: one never can be entirely 
sure whether one has indeed found a state close to the true 
ground state of the system. Hence, to find certifiable methods 
to get ground states seems very timely and important 



FIG. 1: Matrix product states and local hamiltonians. 



plexity of finding ground states using variations over matrix- 
product states as in DMRG. As such, DMRG essentially 
amounts to a local variation of matrices in the matrix-product 
states 10, S] . This is made most explicit in the variant of Ref . 
Jil]. Here, we show that if we allow for a global variation 
over several sites at once - to find the globally best approxi- 
mation in this set and to avoid local optima - one encounters 
a problem which is computationally hard, even in approxima- 
tion. Or, actually more strongly: for any instance of a binary 
quadratic problem (including the NP hard exact satisfiability, 
or the maximum clique or independent set problems) one has 
an identical instance, realized with local translationally invari- 
ant hamiltonians |14]. To prove that this is true, we need 4- 
level systems, 6-local hamiltonians, and a variation over two 
sites. This is much stronger than merely saying that polyno- 
mially constrained problems of high degree as such are com- 
putationally hard, as it could a priori of course well be that 
these hard instances never occur in the specific context un- 
der consideration. Moreover, this reduction can be found with 
polynomial effort. So one faces the ironic and interesting situ- 
ation that when locally varying matrix product states, one has 
an efficient subproblem, but may not get the certifiable true 
ground state of the system. In turn, when aiming at avoiding 
the problem of local minima and varying matrices of several 
sites at once, one encounters a hard problem. 



In this work, we present a first rigorous analysis of the com- 



Local hamiltonians. - The considered hamiltonians are r- 
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local and translationally invariant, 

n— r+l 

up to open boundary conditions (see Fig. 1), where the local 
interaction is governed by some general hamiltonian /i*^*^ = 

T,t,0,...,i=i ha,...,d^^^ ® a^;^'^ ® • • • ® 4*+'^-^)). Here, 
{(Ta''} denotes any local operator basis for site i = 1, . . . , n, 
such as Pauli operators in case of spin-1/2 systems. This in- 
sistence on local hamiltonians renders the assessment of the 
computational complexity fair, as in non-local models, one 
has the freedom to incorporate frustrated higher-dimensional 
systems. 

Matrix product states. - The density matrix renormaliza- 
tion group methods - in several variants - essentially produce 
a sequence of matrix product states (MPS) that better and bet- 
ter approximate the true ground state. These matrix product 
states correspond to non-translationally invariant versions of 
the finitely correlated states [3], which were historically de- 
veloped and studied independently from the DMRG context. 

We consider chains of n sites and d-level constituents, so 
d = 2 for a spin chain. The local basis is denoted by 
{|1), . . . , Then, MPS take the standard form = 

Et...,.=i<^^i:^..^t^|n,^2,...,^n). Here, G 
(C^j x^j+i^ j = 1, . . . , n, are complex matrices. These ma- 
trices, depending on the auxiliary dimension D = maxj Dj 
(the MPS dimension), characterize the MPS. For simplic- 
ity, we impose open boundary conditions. This means that 
D, = Dn+i = 1, i.e., G C^-xi and 

A^^^ , . . . , A^^^ G (D^x^2 are taken to be vectors. We consider 
normalized MPS, meaning that = 1. This is notably 

achieved by using the gauge condition that Yl,i{^^i^)^ — 
1 for each site j llisll . This condition can always be satisfied 
up to similarity transformations JtI], and is in numerical meth- 
ods, most explicitly in Ref. [Isl], insisted upon for numerical 
stability. 

We look at optimizations over MPS, as it is explicitly or 
implicitly being done in DMRG approaches. Needless to say, 
there are many variants. The local variational method for 
finite-size DMRG involves optimizations over single sites Js]], 
essentially equivalent to the B o 5-method of DMRG in case 
of open boundary conditions: For any site j = 1, . . . , n, one 
keeps all matrices A^*^ of the sites i different from j fixed. 
Then, one finds the optimal MPS by varying over the matri- 
ces a[^-' , . . . , A^J-' of site j, satisfying normalization, to mini- 
mize the energy E = {iIj\H\iIj). Then, one takes the next site, 
optimizes over the matrices of that site, and "sweeps", until 
a fixed point is reached 10, 13, Bl . In fact, assertions that in 
gapped systems, DMRG does find the ground state, implicitly 
assume that the globally optimal MPS can be found. 

Main result. - Yet, the problem one actually intends to 
solve is the full problem, so the global optimization problem 
of finding the best matrix product state. Or, if one has an in- 
finite system, one should at least be able to solve the problem 



over several sites at once, in one run, solving the problem over 
these sites, say, of the length scale of the classical correlation 
length. Obviously, we have to allow for sequences of larger 
and larger systems in n, as otherwise, the complexity question 
and the one of finding the ground state no longer makes sense: 
for finite systems, there is always a finite D to exactly write 
out the true ground state. We will now see that there are hard 
instances in the class of variations over several sites at once. 

Problem 1 (DMRG with global variation over several sites) 

Consider a family of translationally invariant hamiltonians 
H of the form as in Eq. ([7]), with n = aD + 6 fa, 6 G R are 
fixed numbers). Let I C {1, . . . , n} some finite subset of sites 
of the hamiltonian. Then, find for any MPS dimension D the 
optimal matrix product state 

1^-)= E A^A^...A^\iui2,...,in) 

by simultaneously varying the matrices A'^^ , • • • , A^J^ of sites 
j G I satisfying the above gauge condition to achieve normal- 
ization (V^IV^) = 1, to minimize the energy E = {ip\H\ip). 

This is the very reasonable and natural variation over several 
sites to minimize the energy, aiming at avoiding local minima. 
Yet, surpringly, one arrives at the subsequent observation: 

Theorem 1 (Hardness of DMRG with global variations) 

Finding the best MPS with varying over several sites is 
NP-hard in D. Moreover, the ground state energy 
can not bel/ D -approximated in polynomial time. 

The latter statement is meant unless P=NP iHIillil]. What 
we aim for is a reduction of this problem to a general binary 
quadratic problem. We start by abstract considerations, and 
then flesh out how we can incorporate this in the MPS set- 
ting. To simplify the notation, we will first consider a setting 
involving N real variables, N even, and will later relate this 
to the dimension D. What follows is not a standard polyno- 
mial reduction to an NP-complete problem: we are heavily 
restricted by the specific form offered by MPS. Yet, it will 
turn out that the freedom that we have is just so sufficient for 
our purposes. We will introduce some new techniques. Read- 
ers only interested in the physical implications may read on 
the with the further discussion of simulatability issues. 

Idea of relating optimization problems. - We start by relat- 
ing a certain class of quadratic continuous problems to a gen- 
eral binary quadratic problem, where variables can only take 
values in {0, 1}, which includes problems as the max clique 
problem. This will be the final form that we want to achieve 
with MPS, when minimizing the energy {ip\H\ip). The ar- 
gument is as follows: For every real {N — 1) x (TV — 1)- 
matrix M with entries having absolute values smaller or equal 
to unity there exists an x TV-matrix Y such that every prob- 
lem of finding the minimal value of hMlF for binary variables 
6i, . . . , ^AT-i G {0, 1} can be written as a problem of find- 
ing the minimal value of xYy^ for xi , . . . , xat , , . . . , ^at G 
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[0, 1]. To see this, we may identify the first — 1 variables 
xi, . . . , xn-1 with 6i, . . . , h^-i- Now, for x, ^ G [0, 1]^, 

N-l N-1 
2 ^ XkVk - ^ (x/e + ^/c) - XnVN 



+ (xi, . . . , Xiv-i)M(xi, . . . , XN-if/{2{N - 1)2) 



takes its minimal value exactly if (i) i/k 



1 



Xk, 



x,y e {0,1}^ (ii) XN 



VN 



1 hold, and (iii) 



(xi, . . . , XAr-i)M(xi, . . . , xn-i)^ takes its minimal value 
over binary variables. Eq. Q can clearly be incorporated in a 
single matrix Y flTll . In this form, we are in the position to ac- 
tually generate exactly this situation in the MPS setting when 
minimizing the energy. Also, the above approximation state- 
ment follows from results in Ref. 1181 , and using estimates for 
the deviation from binary variables, once we are in the posi- 
tion of formulating the problem as a minimization over xYy"^ . 

Incorporating this in MPS. - In order to generate a fair 
worst-case szenario, we are free to take any variational set 
/, fix the other matrices of sites not contained in / appro- 
priately, and may take any local hamiltonian H and local 
physical dimension. Then, energy minimization amounts to 
solving the optimization problem. We will make use of a 
6-local hamiltonian, so h acts non-trivially on 6 subsystems, 
= |i)(i|^6^|2)^2|^6 + |3)(3|^6 + |4)(4|^6foralH.The 
idea now is to think of four "regions" of the matrix product 
chain: To the left and to the right, there will be m = [log2 D] 
sites, forming a "tail" to accumulate the proper range of the 
matrices. The left center consists of 6 sites, and the system 
will be constructed in a way such that the hamiltonian acts 
only non-trivially on these sites. The right center is a chain 
of sites generating "indicator matrices". This means that the 
chain consists for integer N of n = N'^ + 6 + 2m sites, and 
D = 2N'^ -\- N. N labels both the auxiliary matrix dimen- 
sion and the system size (we can always pad the system to get 
n = aD for a, b e N). 

Let us first focus on the left center, embodying 6 sites: 
Sites m + 3 and m + 5 will form the set /, and we keep the 
other matrices fixed, respecting the gauge condition. We take 
^(^+4) ^ ^^iv^ ^^^^ l)/N)eOD-N, and 4""^^^ such that 
the gauge condition is satisfied, the other two matrices of this 
site being zero. Here, we use of the notation E{i^j) for a ma- 
trix, all entries of which are zero, except that E{i^j)ij = 1. 
This matrix has the purpose of selecting appropri- 

ate parts of and ^1^^+^^ ; the matrices and 

^(m+5) ^.jj ^^^^^ incorporate the variables x^ and y^, respec- 
tively. It is not difficult to see that for any complex c/^, dk 
satisfying \ck\ < 1 and \dk\ < I, k = 1,...,7V, we can 
find D X L)-matrices satisfying (^(^+3))t^(^+3) < ^nd 
<li^,suchthat 

(2) 

Conversely, for any solution of Eq. Q we find |c/c|, |(i/c| < 1 
for all k. This follows from exploiting the gauge condi- 



tions for and making use of the specific form of 



(m+3) 



that has been chosen. These numbers ci , . . . , cat and 



(ii, . . . , c^AT are still complex: the key idea we make use of at 
this point is that we can appropriately combine them such that 
only absolute values remain: the role of the binary variables 
will be taken over by Xk = \ck\'^,yk = \dk\'^. 

Generating indicator matrices from matrix products. - To 
the right of the left center, we will append the right center, 
so appropriately chosen N'^ matrices. As we insist on a lo- 
cal hamiltonian, we have no freedom to select only certain 
products: we will always have to deal with all possible prod- 
ucts. We hence have to exploit a certain structure, such that 
from all exponentially many products, we can generate poly- 
nomially many indicator matrices, which are zero except from 
a single non-zero element. These indicator matrices will be 
used in the MPS construction to single out certain elements. 
It is not at all obvious that such matrices generating indicator 
matrices even exist: Yet, for any N x N matrix Y with en- 
tries Yk^i G [0, 1), one can indeed construct D x I)-matrices 
M[^\M!f^ for j = 1, . . . , AT^ with this property For 
every binary word (n, ^2, • • • , ^ats) G {0, 1}^ , we find 



N 



(i) 
+1 



V/, 



k,l 



E(k,l) 




or. 



0, 



for /c, / = 1, . . . , A", where P G R^><^ is defined as P = 
Ylk^i "Z^iLi P^{k + 1, kN + /). Also, each non-zero matrix 
in Eq. (O corresponds to a single binary word. In more col- 
loquial and intuitive terms: we can take an arbitrary product 
of these matrices with lower index 1 or 2 defined by the bi- 
nary word, and multiply it from the left with P, acting here as 
a shift operator. We will then always obtain a matrix with a 
single non-zero element. This single element can be arbitrary, 
forming the desired matrix Y. To check in retrospect that the 
given construction lEoll has this property is straightforward. 

Combining results and minimizing the energy. - We are 
now in the position to put the previous results together, and 
see that the binary optimization problem can indeed be en- 
coded in energy minimization. Concerning the further matri 
ces of the left center, we simply put 

Am+2) 



and 4"^+^^ = Iat © Od-n, A^"^^"' = Oat © 1d-n- 
All other matrices of sites m + 1, m + 2, m + 6 are set to 
zero. The matrices forming the right center are identified 



with A 



(m+6+0 
2+i 



mP, for / 



l,...,Ar2 



1,2, and 



A 



(m+6+O _ _ 



Od foralU = 1,...,A'^ The left 



and right tails are simply padded with matrices such that the 
product of them gives rise to an identity matrix Now, 
we find that we have satisfied the gauge condition for all ma- 
trices of sites {1, . . . , n}\/. Finally, we can combine all this: 
Let the matrices associated with sites {1, . . . , n}\/ be as cho- 
sen above, and the hamiltonian as in Eq. ([T]). Then, energy 
minimization becomes 



(m+3) /i(m+4) . (m+5) 



^1 



^1 



Yk,iE{k,i 



xYy^/N. 



k,l 
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This follows now from a direct evaluation of the overlap of the 
MPS, making use of the above results. This proves the valid- 
ity of the theorem: the optimization problem encountered in 
the variation over matrix-product state is in this form identical 
to solving the respective instance of the binary quadratic prob- 
lem, and for the very same hamiltonian, for every instance of 
this problem, as in an instance of the max clique problem, 
one can find a D such that the variational problem becomes 
identical. This shows that indeed: even within the setting of 
varying over several sites simultaneously when finding opti- 
mal matrix-product states to approximate ground states, one 
encounters computationally difficult NP hard problems. 

Further discussion of the simulatability of quantum many- 
body problems. - In this work, we have addressed the question 
of finding quantum ground states of one-dimensional many- 
body systems on a classical computer, as far as the complex- 
ity of local variations is concerned. This is a question that 
has not explicitly been addressed so far: even if MPS are a 
set faithfully representing the true ground state, how difficult 
is it to find the optimal one? It turns out that in the class 
of problems where one reasonably varies over several sites 
at once contains provably computationally difficult instances, 
even for local one-dimensional hamiltonians. By no means 
should this be read as a statement that DMRG does not work: 
In practice, DMRG obviously gives typically rise to very good 
results. But rather as a warning sign, that to find best approx- 
imations of many-body systems can be computationally hard. 
Moreover, it suggests that to further look for new certifiable 
algorithms to find ground states, including "error bars", at 
least for non-critical systems, should be a very fruitful task. 
A number of questions are implicitly raised here: What is the 
significance of breaking the translational symmetry in MPS 
for ground state approximations? What role does the gauge 
freedom play? Also, what is the exact relationship to QMA 
completeness [23] of local hamiltonian problems? Can good 
bounds be found via polynomial relaxations 11241] ? It would be 
exciting to further look at truly optimal translationally invari- 
ant MPS as approximations of translationally invariant ground 
states, to see under what conditions such approximations are 
truly efficient. It is the hope that the present work further fos- 
ters such considerations. 
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